Fit condition model with environmental and biological predictors

Fit main model (see exploratory scripts and model comparison), visualize results.

Load packages

library(tidyverse); theme_set(theme_light(base_size = 12))
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(viridis)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(gganimate)
library(gifski)
library(latex2exp)
library(patchwork)
library(png)
library(RCurl)
library(wesanderson)
library(qwraps2) 
library(ggcorrplot)
library(ggridges)
# remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)

# To load entire cache in interactive r session, do: 
# qwraps2::lazyload_cache_dir(path = "R/analysis/condition_model_cf_cache/html")

Code for map plots

# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
      axis.text.x = element_text(angle = 90),
      axis.text = element_text(size = 8),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 8),
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
      strip.background = element_rect(fill = "grey95")
      )
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 6),
        strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
        strip.background = element_rect(fill = "gray95"),
        legend.position = c(0.7, 0.02),
        legend.direction = "horizontal"
      )
}

# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2

plot_map_raster <- 
ggplot(swe_coast_proj) + 
  xlim(xmin2, xmax2) +
  ylim(ymin2, ymax2) +
  labs(x = "Longitude", y = "Latitude") +
  geom_sf(size = 0.3) + 
  theme_plot() +
  NULL

plot_map_raster_labels <- 
  plot_map_raster + 
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = 1.9, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = 1.9, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = 1.9, angle = 75)

Read data

# Read the split files and join them
d1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cond_(1_2).csv")
d2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cond_(2_2).csv")

d <- bind_rows(d1, d2)

unique(is.na(d$density_cod))
#> [1] FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE  TRUE

# Calculate standardized variables
d <- d %>% 
  mutate(ln_length_cm = log(length_cm),
         ln_weight_g = log(weight_g),
         year_ct = year - mean(year),
         biomass_her_sc = biomass_her,
         biomass_her_sd_sc = biomass_her_sd,
         biomass_spr_sc = biomass_spr,
         biomass_spr_sd_sc = biomass_spr_sd,
         density_cod_sc = density_cod,
         density_cod_rec_sc = density_cod_rec,
         density_fle_sc = density_fle,
         density_fle_rec_sc = density_fle_rec,
         density_saduria_sc = density_saduria,
         density_saduria_rec_sc = density_saduria_rec,
         depth_sc = depth,
         depth_rec_sc = depth_rec,
         depth_sd_sc = depth_sd,
         oxy_sc = oxy,
         oxy_rec_sc = oxy_rec,
         oxy_sd_sc = oxy_sd,
         temp_sc = temp,
         temp_rec_sc = temp_rec,
         temp_sd_sc = temp_sd,
         year_f = as.factor(year))  %>%
  mutate_at(c("biomass_her_sc", "biomass_her_sd_sc", "biomass_spr_sc", "biomass_spr_sd_sc",
              "density_cod_sc", "density_cod_rec_sc", 
              "density_fle_sc", "density_fle_rec_sc", 
              "density_saduria_sc", "density_saduria_rec_sc", 
              "depth_sc", "depth_rec_sc", "depth_sd_sc",
              "oxy_sc", "oxy_rec_sc", "oxy_sd_sc",
              "temp_sc", "temp_rec_sc", "temp_sd_sc"
              ),
            ~(scale(.) %>% as.vector)) %>% 
  mutate(year = as.integer(year))

unique(is.na(d))
#>      weight_g length_cm  year   lat   lon ices_rect sub_div depth   oxy  temp
#> [1,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [2,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [3,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [4,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [5,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [6,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#> [7,]    FALSE     FALSE FALSE FALSE FALSE     FALSE   FALSE FALSE FALSE FALSE
#>          X     Y density_cod density_fle density_cod_data density_fle_data
#> [1,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [2,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [3,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [4,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#> [5,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#> [6,] FALSE FALSE       FALSE       FALSE            FALSE            FALSE
#> [7,] FALSE FALSE       FALSE       FALSE             TRUE             TRUE
#>      density_saduria density_saduria_rec density_cod_rec density_fle_rec
#> [1,]           FALSE               FALSE           FALSE           FALSE
#> [2,]           FALSE               FALSE           FALSE           FALSE
#> [3,]           FALSE               FALSE           FALSE           FALSE
#> [4,]           FALSE               FALSE           FALSE           FALSE
#> [5,]           FALSE                TRUE            TRUE            TRUE
#> [6,]           FALSE                TRUE            TRUE            TRUE
#> [7,]           FALSE               FALSE           FALSE           FALSE
#>      depth_rec oxy_rec temp_rec biomass_spr biomass_her density_saduria_sd
#> [1,]     FALSE   FALSE    FALSE       FALSE       FALSE              FALSE
#> [2,]     FALSE   FALSE    FALSE        TRUE        TRUE              FALSE
#> [3,]     FALSE   FALSE    FALSE        TRUE        TRUE              FALSE
#> [4,]     FALSE   FALSE    FALSE       FALSE       FALSE              FALSE
#> [5,]      TRUE    TRUE     TRUE        TRUE        TRUE              FALSE
#> [6,]      TRUE    TRUE     TRUE        TRUE        TRUE              FALSE
#> [7,]     FALSE   FALSE    FALSE        TRUE        TRUE              FALSE
#>      density_cod_sd density_fle_sd depth_sd oxy_sd temp_sd biomass_spr_sd
#> [1,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [2,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [3,]          FALSE          FALSE    FALSE  FALSE   FALSE           TRUE
#> [4,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [5,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [6,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#> [7,]          FALSE          FALSE    FALSE  FALSE   FALSE          FALSE
#>      biomass_her_sd ln_length_cm ln_weight_g year_ct biomass_her_sc
#> [1,]          FALSE        FALSE       FALSE   FALSE          FALSE
#> [2,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [3,]           TRUE        FALSE       FALSE   FALSE           TRUE
#> [4,]          FALSE        FALSE       FALSE   FALSE          FALSE
#> [5,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [6,]          FALSE        FALSE       FALSE   FALSE           TRUE
#> [7,]          FALSE        FALSE       FALSE   FALSE           TRUE
#>      biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc density_cod_sc
#> [1,]             FALSE          FALSE             FALSE          FALSE
#> [2,]             FALSE           TRUE             FALSE          FALSE
#> [3,]              TRUE           TRUE              TRUE          FALSE
#> [4,]             FALSE          FALSE             FALSE          FALSE
#> [5,]             FALSE           TRUE             FALSE          FALSE
#> [6,]             FALSE           TRUE             FALSE          FALSE
#> [7,]             FALSE           TRUE             FALSE          FALSE
#>      density_cod_rec_sc density_fle_sc density_fle_rec_sc density_saduria_sc
#> [1,]              FALSE          FALSE              FALSE              FALSE
#> [2,]              FALSE          FALSE              FALSE              FALSE
#> [3,]              FALSE          FALSE              FALSE              FALSE
#> [4,]              FALSE          FALSE              FALSE              FALSE
#> [5,]               TRUE          FALSE               TRUE              FALSE
#> [6,]               TRUE          FALSE               TRUE              FALSE
#> [7,]              FALSE          FALSE              FALSE              FALSE
#>      density_saduria_rec_sc depth_sc depth_rec_sc depth_sd_sc oxy_sc oxy_rec_sc
#> [1,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#> [2,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#> [3,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#> [4,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#> [5,]                   TRUE    FALSE         TRUE       FALSE  FALSE       TRUE
#> [6,]                   TRUE    FALSE         TRUE       FALSE  FALSE       TRUE
#> [7,]                  FALSE    FALSE        FALSE       FALSE  FALSE      FALSE
#>      oxy_sd_sc temp_sc temp_rec_sc temp_sd_sc year_f
#> [1,]     FALSE   FALSE       FALSE      FALSE  FALSE
#> [2,]     FALSE   FALSE       FALSE      FALSE  FALSE
#> [3,]     FALSE   FALSE       FALSE      FALSE  FALSE
#> [4,]     FALSE   FALSE       FALSE      FALSE  FALSE
#> [5,]     FALSE   FALSE        TRUE      FALSE  FALSE
#> [6,]     FALSE   FALSE        TRUE      FALSE  FALSE
#> [7,]     FALSE   FALSE       FALSE      FALSE  FALSE

unique(is.na(d$density_cod_rec))
#> [1] FALSE  TRUE
unique(is.na(d$density_cod))
#> [1] FALSE

# Drop NA covariates for the correlation plot
d <- d %>% drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
                   density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc, 
                   density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
                   oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc)



# Plot correlation between variables
d_cor <- d %>% dplyr::select("biomass_her_sc", "biomass_her_sd_sc",
                             "biomass_spr_sc", "biomass_spr_sd_sc",
                             "density_cod_sc", "density_cod_rec_sc", 
                             "density_fle_sc", "density_fle_rec_sc", 
                             "density_saduria_sc", "density_saduria_rec_sc", 
                             "depth_sc", "depth_rec_sc",
                             "oxy_sc", "oxy_rec_sc",
                             "temp_sc", "temp_rec_sc")

corr <- round(cor(d_cor), 1)

ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.3))

ggsave("figures/supp/condition/correlation_variables.png", width = 6.5, height = 6.5, dpi = 600)

Read the prediction grids

pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")

pred_grid <- bind_rows(pred_grid1, pred_grid2)

unique(is.na(pred_grid$density_cod))
#> [1] FALSE
unique(is.na(pred_grid$density_cod_rec))
#> [1] FALSE

pred_grid <- pred_grid %>% mutate(year = as.integer(year),
                                  year_f = as.factor(year))

# Scale the variables with respect to data! First calculate means in data
data_means <- d %>%
  dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
                density_cod, density_cod_rec, 
                density_fle, density_fle_rec, 
                density_saduria, density_saduria_rec, 
                depth, depth_rec, depth_sd,
                oxy, oxy_rec, oxy_sd,
                temp, temp_rec, temp_sd) %>%
  mutate_all(~mean(.)) %>%
  distinct(.keep_all = TRUE)

data_stdev <- d %>%
  dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
                density_cod, density_cod_rec, 
                density_fle, density_fle_rec, 
                density_saduria, density_saduria_rec, 
                depth, depth_rec, 
                oxy, oxy_rec, 
                temp, temp_rec
                ) %>%
  mutate_all(~sd(.)) %>%
  distinct(.keep_all = TRUE)

# Before actually scaling, replace the ices-rectangle pelagic values that are NA with the mean across rectangles 
# in the sub division. Replace the sub-division pelagic values that are NA with the mean in the year.
# Note that the sub-division values are not the sum of the rectangles (due to missing rectangles), so 
# I need to calculate a sub-division mean across rectangles within each sub-division
pred_grid <- pred_grid %>% 
  group_by(year) %>% 
  mutate(median_biomass_sprat_across_sub_div = median(biomass_spr_sd, na.rm = TRUE),
         median_biomass_herring_across_sub_div = median(biomass_her_sd, na.rm = TRUE)) %>% 
  ungroup() %>% # Replace sub_divsion NA's with the median across sub_divisions in that year
  mutate(biomass_spr_sd = ifelse(is.na(biomass_spr_sd) == TRUE, median_biomass_sprat_across_sub_div, biomass_spr_sd),
         biomass_her_sd = ifelse(is.na(biomass_her_sd) == TRUE, median_biomass_herring_across_sub_div, biomass_her_sd)) %>% 
  group_by(year, sub_div) %>% 
  mutate(median_biomass_sprat_across_rect_in_sub_div = median(biomass_spr, na.rm = TRUE),
         median_biomass_herring_across_rect_in_sub_div = median(biomass_her, na.rm = TRUE)) %>% 
  ungroup() %>% # Replace rectangle NA's with the median across rectangles in that year and sub-division
  mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, median_biomass_sprat_across_rect_in_sub_div, biomass_spr),
         biomass_her = ifelse(is.na(biomass_her) == TRUE, median_biomass_herring_across_rect_in_sub_div, biomass_her)) %>% 
  # Since I still have some NAs (some sub-divisions do not have a single rectangle in some years), I will will those rectangles 
  # with the sub-division value divided by the number of rectangles in that sub division
  group_by(year, sub_div) %>% 
  mutate(n_rect = length(unique(ices_rect))) %>% 
  ungroup() %>% 
  mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, biomass_spr_sd/n_rect, biomass_spr),
         biomass_her = ifelse(is.na(biomass_her) == TRUE, biomass_her_sd/n_rect, biomass_her))

pred_grid <- pred_grid %>%
  mutate(ln_length_cm = log(1),
         year_ct = year - mean(year),
         biomass_her_sc = (biomass_her - data_means$biomass_her) / data_stdev$biomass_her,
         biomass_her_sd_sc = (biomass_her_sd - data_means$biomass_her_sd) / data_stdev$biomass_her_sd,
         biomass_spr_sc = (biomass_spr - data_means$biomass_spr) / data_stdev$biomass_spr,
         biomass_spr_sd_sc = (biomass_spr_sd - data_means$biomass_spr_sd) / data_stdev$biomass_spr_sd,
         density_cod_sc = (density_cod - data_means$density_cod) / data_stdev$density_cod,
         density_cod_rec_sc = (density_cod_rec - data_means$density_cod_rec) / data_stdev$density_cod_rec,
         density_fle_sc = (density_fle - data_means$density_fle) / data_stdev$density_fle,
         density_fle_rec_sc = (density_fle_rec - data_means$density_fle_rec) / data_stdev$density_fle_rec,
         density_saduria_sc = (density_saduria - data_means$density_saduria) / data_stdev$density_saduria,
         density_saduria_rec_sc = (density_saduria_rec - data_means$density_saduria_rec) / data_stdev$density_saduria_rec,
         depth_sc = (depth - data_means$depth) / data_stdev$depth,
         depth_rec_sc = (depth_rec - data_means$depth_rec) / data_stdev$depth_rec,
         oxy_sc = (oxy - data_means$oxy) / data_stdev$oxy,
         oxy_rec_sc = (oxy_rec - data_means$oxy_rec) / data_stdev$oxy_rec,
         temp_sc = (temp - data_means$temp) / data_stdev$temp,
         temp_rec_sc = (temp_rec - data_means$temp_rec) / data_stdev$temp_rec,
         )

# Plot on map:
ggplot(pred_grid2, aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

pred_grid %>% 
  drop_na() %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

pred_grid %>% 
  drop_na(biomass_spr_sd) %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

unique(is.na(pred_grid))
#>          X     Y depth  year   oxy  temp   lon   lat sub_div ices_rect
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE   FALSE     FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE   FALSE     FALSE
#>      density_saduria biomass_spr biomass_her biomass_spr_sd biomass_her_sd
#> [1,]           FALSE       FALSE       FALSE          FALSE          FALSE
#> [2,]           FALSE       FALSE       FALSE          FALSE          FALSE
#>      density_cod density_fle depth_rec temp_rec oxy_rec density_cod_rec
#> [1,]       FALSE       FALSE     FALSE    FALSE   FALSE           FALSE
#> [2,]       FALSE       FALSE     FALSE    FALSE   FALSE           FALSE
#>      density_fle_rec density_saduria_rec depth_sd temp_sd oxy_sd density_cod_sd
#> [1,]           FALSE               FALSE    FALSE   FALSE  FALSE          FALSE
#> [2,]           FALSE               FALSE    FALSE   FALSE  FALSE          FALSE
#>      density_fle_sd density_saduria_sd year_f
#> [1,]          FALSE              FALSE  FALSE
#> [2,]          FALSE              FALSE  FALSE
#>      median_biomass_sprat_across_sub_div median_biomass_herring_across_sub_div
#> [1,]                               FALSE                                 FALSE
#> [2,]                               FALSE                                 FALSE
#>      median_biomass_sprat_across_rect_in_sub_div
#> [1,]                                       FALSE
#> [2,]                                        TRUE
#>      median_biomass_herring_across_rect_in_sub_div n_rect ln_length_cm year_ct
#> [1,]                                         FALSE  FALSE        FALSE   FALSE
#> [2,]                                          TRUE  FALSE        FALSE   FALSE
#>      biomass_her_sc biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc
#> [1,]          FALSE             FALSE          FALSE             FALSE
#> [2,]          FALSE             FALSE          FALSE             FALSE
#>      density_cod_sc density_cod_rec_sc density_fle_sc density_fle_rec_sc
#> [1,]          FALSE              FALSE          FALSE              FALSE
#> [2,]          FALSE              FALSE          FALSE              FALSE
#>      density_saduria_sc density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc
#> [1,]              FALSE                  FALSE    FALSE        FALSE  FALSE
#> [2,]              FALSE                  FALSE    FALSE        FALSE  FALSE
#>      oxy_rec_sc temp_sc temp_rec_sc
#> [1,]      FALSE   FALSE       FALSE
#> [2,]      FALSE   FALSE       FALSE

pred_grid %>% 
  drop_na(biomass_spr_sc, biomass_spr_sd_sc, biomass_her_sc, biomass_her_sd_sc) %>% 
  ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)

Make spde mesh

spde <- make_mesh(d, xy_cols = c("X", "Y"),
                  n_knots = 100,  
                  type = "kmeans", seed = 42)

# Plot and save spde
png(file = "figures/supp/condition/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()

Fit models

Non spatial model to get the average length weight relationship. The main model’s response variable is the ratio of the observed vs this predicted average weight.

ptm <- proc.time()
mnull <- sdmTMB(
  formula = ln_weight_g ~ ln_length_cm,
  data = d,
  time = NULL,
  mesh = spde,
  family = student(link = "identity", df = 5),
  spatiotemporal = "off",
  spatial = "off",
  silent = TRUE,
  reml = FALSE
  )
proc.time() - ptm
#>    user  system elapsed 
#>   1.945   0.263   2.765

# Extract average a and b across the total dataset
a <- tidy(mnull)[1, "estimate"]
b <- tidy(mnull)[2, "estimate"]
a
#> [1] -4.638063
b
#> [1] 2.983866
d$weight_g_avg_pred <- exp(a)*d$length_cm^b

plot(d$weight_g ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$weight_g); abline(a = 0, b = 1, col = "red")

d$le_cren <- d$weight_g / d$weight_g_avg_pred
d$log_le_cren <- log(d$le_cren)
hist(d$le_cren)

Compare breakpoint vs linear oxygen effect in simplified model

mnull_oxy_bp <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + breakpt(oxy_sc) + oxy_rec_sc,
  time_varying = NULL,
  data = d,
  time = "year",
  mesh = spde,
  family = student(link = "identity", df = 5), 
  spatiotemporal = "iid",
  spatial = "on",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1))

mnull_oxy <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + oxy_sc + oxy_rec_sc,
  time_varying = NULL,
  data = d,
  time = "year",
  mesh = spde,
  family = student(link = "identity", df = 5), 
  spatiotemporal = "iid",
  spatial = "on",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1))

AIC(mnull_oxy_bp, mnull_oxy)
min(AIC(mnull_oxy_bp), AIC(mnull_oxy)) 
#> [1] -144518.8

sanity(mnull_oxy_bp)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
sanity(mnull_oxy)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large

tidy(mnull_oxy)
tidy(mnull_oxy_bp)

nd_bp_simple <- data.frame(
  oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100), 
  year = 2019L,
  year_f = as.factor(2019),
  oxy_rec_sc = 0
)

p_bp_simple <- predict(mnull_oxy_bp, newdata = nd_bp_simple,
                       se_fit = TRUE, re_form = NA, xy_cols = c("X", "Y"))

ggplot(p_bp_simple, aes(oxy_sc, est, 
  ymin = est - 1.96 * est_se, 
  ymax = est + 1.96 * est_se)) +
  geom_line() + geom_ribbon(alpha = 0.4)

ggplot(d, aes(oxy_sc, oxy)) + 
  geom_point() + 
  geom_hline(yintercept = (1.0700*sd(d$oxy)) + mean(d$oxy)) +
  geom_vline(xintercept = 1.0700)

Fit main model

ptm <- proc.time()
mfull <- sdmTMB(
  formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
    biomass_spr_sc + biomass_spr_sd_sc +
    density_cod_sc + density_cod_rec_sc + 
    density_fle_sc + density_fle_rec_sc + 
    density_saduria_sc + density_saduria_rec_sc + 
    depth_sc + depth_rec_sc + 
    oxy_sc + oxy_rec_sc + 
    temp_sc + temp_rec_sc,
    data = d,
  time = "year",
  mesh = spde, 
  family = student(link = "identity", df = 5),
  spatiotemporal = "iid",
  spatial = "on",
  silent = TRUE,
  reml = FALSE,
  control = sdmTMBcontrol(newton_loops = 1)
  )
proc.time() - ptm
#>    user  system elapsed 
#> 721.928  20.345 747.627

tidy(mfull, conf.int = TRUE) %>% arrange(desc(abs(estimate)))
tidy(mfull, effects = "ran_par", conf.int = TRUE) 
# Check model convergence
sanity(mfull)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
mcmc_res <- residuals(mfull, type = "mle-mcmc", mcmc_iter = 5000, mcmc_warmup = 2000)

png(file = "figures/supp/condition/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()

d$residuals <- mcmc_res[, 1]

Calculate approximate r^2 manually using a modifed r2.sdmTMB

b_fe <- tidy(mfull)
b_fe <- stats::setNames(b_fe$estimate, b_fe$term)
X <- mfull$tmb_data$X_ij
X <- matrix(unlist(X), ncol = length(b_fe))

VarF <- var(as.vector(b_fe %*% t(X))) # variance from fixed-effects
 
b_ran <- tidy(mfull, "ran_par")
sigma_O <- b_ran$estimate[b_ran$term == "sigma_O"] # spatial standard deviation
sigma_E <- b_ran$estimate[b_ran$term == "sigma_E"] # spatiotemporal standard deviation
VarO <- sigma_O^2 # spatial variance
VarE <- sigma_E^2 # spatiotemporal variance
 
# Calculate obs - pred
d_r2 <- d
d_r2$pred <- predict(mfull)
d_r2$resid_manual <- d_r2$le_cren - d_r2$pred$est
#hist(d_r2$resid_manual)

VarR <- var(as.vector(d_r2$resid_manual)) # "residual" variance

# Marginal R2s
denominator <- VarF + VarO + VarE + VarR

marg <- VarF/denominator

out <- data.frame(
    marginal = marg,
    partial_spatial = VarO/denominator,
    partial_spatiotemporal = VarE/denominator,
    conditional = (VarF + VarO + VarE)/denominator,
    all_random = (VarO + VarE)/denominator,
    marginal_random_ratio = marg / ((VarO + VarE)/denominator),
    random_marginal_ratio = ((VarO + VarE)/denominator) / marg
  )

out

# What if I skip the year effect?
b_fe2 <- tidy(mfull) %>% filter(!grepl('year', term))
#> filter: removed 27 rows (63%), 16 rows remaining
b_fe2 <- stats::setNames(b_fe2$estimate, b_fe2$term)

X2 <- mfull$tmb_data$X_ij
X2 <- matrix(unlist(X2), ncol = length(b_fe))[, 28:43] # skip year columns
VarF2 <- var(as.vector(b_fe2 %*% t(X2))) # variance from fixed-effects

denominator3 <- VarF2 + VarO + VarE + VarR

marg3 <- VarF2/denominator3
marg3
#> [1] 0.05794511

Plot residuals

# Residuals on map
plot_map_raster +
  geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals)) +
  scale_colour_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  labs(color = "residuals") +
  NULL


ggsave("figures/supp/condition/spatial_resid.png", width = 6.5, height = 6.5, dpi = 600)

# Residuals vs length and year
ggplot(d, aes(x = length_cm, y = residuals, color = residuals)) +
  geom_point(size = 0.1, alpha = 1) +
  geom_hline(yintercept = 0, size = 0.5, color = "black", linetype = 2, alpha = 0.5) +
  labs(x = "length [cm]", y = "Residuals", color = "") +
  stat_smooth(color = "black", size = 0.5) +
  facet_wrap(~year, ncol = 5) +
  scale_color_gradient2() +
  theme_facet_map()


ggsave("figures/supp/condition/residuals_vs_length_and_year.png", width = 6.5, height = 6.5, dpi = 600)

ggplot(d, aes(year, length_cm, z = residuals)) +
  stat_summary_2d(bins = 40) +
  scale_fill_gradient2()


# ggsave("figures/supp/condition/residuals_vs_length_and_year2.png", width = 6.5, height = 6.5, dpi = 600)

Annual index using get_index_sims method to avoid slow bias correction

ncells <- filter(pred_grid, year == max(pred_grid$year)) %>% nrow()
#> filter: removed 240,214 rows (96%), 9,239 rows remaining

pred_avg_sim <- predict(mfull, newdata = pred_grid, nsim = 100)

index_avg_sim <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_sims <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)),
                                 return_sims = TRUE) # This is for calculation on samples!
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

write.csv(index_avg_sim, "output/index_avg_sim.csv")

Evaluate sensitivity to using a pred grid with a subset excluding low density areas

# Now do the same but excluding the areas not sampled in the pred grid
ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 20 & depth < 120) %>% nrow()
#> filter: removed 242,779 rows (97%), 6,674 rows remaining
pred_grid_sub <- filter(pred_grid, depth > 20 & depth < 120)
#> filter: removed 69,255 rows (28%), 180,198 rows remaining

pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 100)

index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(1/ncells_sub, nrow(pred_avg_sim_sub)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

# Combine and plot & compare
index_avg_sim_comp <- bind_rows(mutate(index_avg_sim, area = "full"),
                                mutate(index_avg_sim_sub, area = "subset"))
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'area' (character) with one unique value and 0% NA

ggplot(index_avg_sim_comp, aes(year, log(est), ymin = log(lwr), ymax = log(upr), color = area)) +
  geom_line() +
  geom_ribbon(alpha = 0.4)

ggplot(index_avg_sim_comp, aes(year, est, ymin = lwr, ymax = upr, color = area)) +
  geom_line() +
  geom_ribbon(alpha = 0.4)

Difference in Le Cren condition factor between the end and start of time period using sims

colnames(pred_grid)
#>  [1] "X"                                            
#>  [2] "Y"                                            
#>  [3] "depth"                                        
#>  [4] "year"                                         
#>  [5] "oxy"                                          
#>  [6] "temp"                                         
#>  [7] "lon"                                          
#>  [8] "lat"                                          
#>  [9] "sub_div"                                      
#> [10] "ices_rect"                                    
#> [11] "density_saduria"                              
#> [12] "biomass_spr"                                  
#> [13] "biomass_her"                                  
#> [14] "biomass_spr_sd"                               
#> [15] "biomass_her_sd"                               
#> [16] "density_cod"                                  
#> [17] "density_fle"                                  
#> [18] "depth_rec"                                    
#> [19] "temp_rec"                                     
#> [20] "oxy_rec"                                      
#> [21] "density_cod_rec"                              
#> [22] "density_fle_rec"                              
#> [23] "density_saduria_rec"                          
#> [24] "depth_sd"                                     
#> [25] "temp_sd"                                      
#> [26] "oxy_sd"                                       
#> [27] "density_cod_sd"                               
#> [28] "density_fle_sd"                               
#> [29] "density_saduria_sd"                           
#> [30] "year_f"                                       
#> [31] "median_biomass_sprat_across_sub_div"          
#> [32] "median_biomass_herring_across_sub_div"        
#> [33] "median_biomass_sprat_across_rect_in_sub_div"  
#> [34] "median_biomass_herring_across_rect_in_sub_div"
#> [35] "n_rect"                                       
#> [36] "ln_length_cm"                                 
#> [37] "year_ct"                                      
#> [38] "biomass_her_sc"                               
#> [39] "biomass_her_sd_sc"                            
#> [40] "biomass_spr_sc"                               
#> [41] "biomass_spr_sd_sc"                            
#> [42] "density_cod_sc"                               
#> [43] "density_cod_rec_sc"                           
#> [44] "density_fle_sc"                               
#> [45] "density_fle_rec_sc"                           
#> [46] "density_saduria_sc"                           
#> [47] "density_saduria_rec_sc"                       
#> [48] "depth_sc"                                     
#> [49] "depth_rec_sc"                                 
#> [50] "oxy_sc"                                       
#> [51] "oxy_rec_sc"                                   
#> [52] "temp_sc"                                      
#> [53] "temp_rec_sc"

# Plot the sims
ggplot(index_avg_sims, aes(y = as.factor(year), x = .value)) + 
  geom_density_ridges_gradient(rel_min_height = 0.01, alpha = 0.8, color = "grey10", size = 0.2,
                               quantile_lines = TRUE, quantiles = 0.5) + 
  coord_flip()
#> Picking joint bandwidth of 0.00302


# index_avg_sims %>% 
#   filter(year %in% c(1993, 2019)) %>% 
#   ggplot(., aes(as.factor(year), .value)) +
#   geom_violin(width = 1.3)

# Now reshape data again to calculate row-wise differences
year_diff <- index_avg_sims %>% 
  filter(year %in% c(1993, 2019)) %>% 
  #mutate(.value = log(.value)) %>% 
  pivot_wider(names_from = year, values_from = .value) %>% 
  mutate("2019-1993" = `2019` - `1993`) %>% 
  pivot_longer(2:4) %>% 
  ungroup()
#> filter: removed 2,500 rows (93%), 200 rows remaining
#> pivot_wider: reorganized (year, .value) into (1993, 2019) [was 200x3, now 100x3]
#> mutate: new variable '2019-1993' (double) with 100 unique values and 0% NA
#> pivot_longer: reorganized (1993, 2019, 2019-1993) into (name, value) [was 100x4, now 300x3]
#> ungroup: no grouping variables

# Plot these
p1 <- ggplot(year_diff) +
  geom_density(data = filter(year_diff, !name == "2019-1993"),
               aes(value, fill = name), alpha = 0.6) +
  scale_fill_brewer(palette = "Dark2") +
  coord_cartesian(expand = 0)
#> filter: removed 100 rows (33%), 200 rows remaining

p2 <- ggplot(year_diff) +
  geom_density(data = filter(year_diff, name == "2019-1993"),
               aes(value), alpha = 0.6, fill = "grey") +
  geom_vline(xintercept = 0, linetype = 2) 
#> filter: removed 200 rows (67%), 100 rows remaining

p1/p2


# Calculate quantiles for each level (to go in the results section)
year_diff %>% 
  group_by(name) %>% 
  summarise(median = quantile(value, probs = 0.5),
            upr97.5 = quantile(value, probs = 0.975),
            lwr2.5 = quantile(value, probs = 0.025))
#> group_by: one grouping variable (name)
#> summarise: now 3 rows and 4 columns, ungrouped

# Percent change in condition factor
year_diff %>% 
  group_by(name) %>% 
  pivot_wider(names_from = name, values_from = value) %>%
  mutate(perc_change = ((`2019`-`1993`)/`1993`)*100) %>% 
  ungroup() %>% 
  summarise(lwr2.5 = quantile(perc_change, probs = 0.025),
            median = quantile(perc_change, probs = 0.5),
            upr97.5 = quantile(perc_change, probs = 0.975))
#> group_by: one grouping variable (name)
#> pivot_wider: reorganized (name, value) into (1993, 2019, 2019-1993) [was 300x3, now 100x4]
#> mutate: new variable 'perc_change' (double) with 100 unique values and 0% NA
#> ungroup: no grouping variables
#> summarise: now one row and 3 columns, ungrouped

Simulated annual condition factor by sub-division

# From these models, predict annual condition factor
# I will use the prediction grid with ALL covariates set to 0, incl. depth, temp and oxygen
# Grabbing the number of cells to help with calculating the average
ncells24 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 24) %>% nrow()
ncells25 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 25) %>% nrow()
ncells26 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 26) %>% nrow()
ncells27 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 27) %>% nrow()
ncells28 <- filter(pred_grid, year == max(pred_grid$year) & sub_div == 28) %>% nrow()

# Use the `area` argument here to turn the total into an average by giving it one over the number of cells
pred_avg_24_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 24), sims = 100)
#> Warning: The `sims` argument of `predict.sdmTMB()` is deprecated as of sdmTMB 0.0.21.
#> Please use the `nsim` argument instead.
pred_avg_25_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 25), sims = 100)
pred_avg_26_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 26), sims = 100)
pred_avg_27_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 27), sims = 100)
pred_avg_28_sim <- predict(mfull, newdata = filter(pred_grid, sub_div == 28), sims = 100)

index_avg_24_sim <- get_index_sims(pred_avg_24_sim, area = rep(1/ncells24, nrow(pred_avg_24_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_25_sim <- get_index_sims(pred_avg_25_sim, area = rep(1/ncells25, nrow(pred_avg_25_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_26_sim <- get_index_sims(pred_avg_26_sim, area = rep(1/ncells26, nrow(pred_avg_26_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_27_sim <- get_index_sims(pred_avg_27_sim, area = rep(1/ncells27, nrow(pred_avg_27_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
index_avg_28_sim <- get_index_sims(pred_avg_28_sim, area = rep(1/ncells28, nrow(pred_avg_28_sim)))
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.

# Combine and plot
index_avg_24_sim <- index_avg_24_sim %>% mutate(sub_div = 24)
index_avg_25_sim <- index_avg_25_sim %>% mutate(sub_div = 25)
index_avg_26_sim <- index_avg_26_sim %>% mutate(sub_div = 26)
index_avg_27_sim <- index_avg_27_sim %>% mutate(sub_div = 27)
index_avg_28_sim <- index_avg_28_sim %>% mutate(sub_div = 28)

div_index <- bind_rows(index_avg_24_sim, index_avg_25_sim, index_avg_26_sim, index_avg_27_sim, index_avg_28_sim)

Condition factor over time and individual distributions

div_index2 <- bind_rows(mutate(div_index, sub_div = as.factor(sub_div)),
                        mutate(index_avg_sim, sub_div = as.factor("Total")))

pal <- c(brewer.pal(n = 5, name = "Dark2"), "grey10")

# Remember to log now that we use sims
cond_index <- ggplot(div_index2, aes(year, est, color = sub_div, linetype = sub_div)) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  labs(y = "Predicted Le Cren's condition factor", x = "Year") +
  geom_line(alpha = 0.8) +
  geom_point(alpha = 0) + 
  scale_color_manual(values = pal, name = "ICES subdivision") +
  scale_linetype_manual(values = c(1,1,1,1,1,0)) +
  geom_point(data = index_avg_sim, aes(year, est), inherit.aes = FALSE, 
             size = 2.5, shape = 21, fill = "gray1", color = "white") +
  geom_errorbar(data = index_avg_sim, inherit.aes = FALSE, 
                aes(x = year, ymax = upr, ymin = lwr),
                width = 0, alpha = 0.8, color = "gray1", size = 0.7) +
  theme(axis.title = element_text(size = 9),
        legend.position = "bottom",
        legend.background = element_blank(),
        legend.title = element_text(size = 9),
        legend.spacing.x = unit(0, 'cm'),
        plot.margin = unit(c(0.4, 0.4, 0.4, 0), "cm"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  guides(linetype = "none",
         color = guide_legend(title.position = "top", title.hjust = 0.5, nrow = 1,
                              keywidth = 0.1, keyheight = 0.1, default.unit = "inch",
                              override.aes = list(linetype = c(1,1,1,1,1,0),
                                                  shape = c(NA,NA,NA,NA,NA,16),
                                                  alpha = rep(0.8, 6),
                                                  color = pal))) + 
  NULL

cond_index


# Ridge-plot of data
data_ridge <- ggplot(d, aes(x = le_cren, y = fct_rev(as.factor(year)), fill = factor(stat(quantile)))) +
  stat_density_ridges(
    rel_min_height = .01,
    scale = 3,
    geom = "density_ridges_gradient",
    calc_ecdf = TRUE,
    alpha = 0.1,
    quantiles = c(0.25, 0.75)) +
  scale_fill_manual(
    name = "Probability",
    values = c("grey50", "gray70", "grey90"),
    labels = c("[0,0.25]", "[0.25,0.75]", "[0.75,1]")) +
  xlim(0.65, 1.55) + 
  labs(x = "Le Cren's condidition factor",
       y = "Year") + 
  guides(fill = guide_legend(nrow = 1, title.position = "top", title.hjust = 0.5)) + 
  theme(axis.title = element_text(size = 9),
        legend.position = "bottom",
        legend.key.size = unit(0.5, "line"),
        legend.background = element_blank(),
        legend.title = element_text(size = 9),
        legend.spacing.x = unit(0.1, 'cm'),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        plot.margin = unit(c(0.4, 0.8, 0.4, 0), "cm")) +
  NULL

data_ridge + cond_index + plot_annotation(tag_levels = 'A') +
   plot_layout(widths = c(1, 2.2))
#> Warning: Removed 926 rows containing non-finite values (stat_density_ridges).


ggsave("figures/cond_index.pdf", width = 17, height = 17, units = "cm")
#> Warning: Removed 926 rows containing non-finite values (stat_density_ridges).


## Calculate delta_fulton and relate to delta_oxygen, by subdivision
div_index_delta_cond <- div_index2 %>% 
  group_by(sub_div) %>% 
  filter(year %in% c(1993, 1994, 1995, 2017, 2018, 2019)) %>% 
  dplyr::select(year, log_est, sub_div) %>% 
  pivot_wider(names_from = year, values_from = log_est) %>% 
  mutate(start_cond = mean(`1993`, `1994`, `1995`),
         end_cond = mean(`2017`, `2018`, `2019`),
         delta_cond = end_cond - start_cond) %>% 
  ungroup() %>% 
  dplyr::select(sub_div, delta_cond)

write.csv(div_index_delta_cond, "output/div_index_delta_cond.csv")

Effect sizes

# What is one standard deviation of oxygen?
sd(d$oxy)
#> [1] 1.846734

# Extract random and fixed coefficients from the full model
mfull_est <- bind_rows(tidy(mfull, effects = "ran_par", conf.int = TRUE) %>%
                         filter(term %in% c("sigma_O", "sigma_E")),
                       
                       tidy(mfull, effects = "fixed", conf.int = TRUE)  %>% 
                         filter(!grepl('year', term))) %>%
  
  mutate(term = factor(term))

mfull_est <- mfull_est %>% 
  mutate(group = "Herring",
         group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
         group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
         group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
         group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
         group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
         group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
         group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
         group = ifelse(term == "sigma_E", "Random", group),
         group = ifelse(term == "sigma_O", "Random", group)
         )

mfull_est <- mfull_est %>% 
  mutate(scale = "Sub-division",
         scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
                                    "density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
         scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
                                    "depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
         scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
         scale = ifelse(term == "sigma_O", "Spatial/spatiotemporal s.d.", scale)
         )

# Plot effects
# This is the order changed labels should go in!
levels(mfull_est$term)
#>  [1] "biomass_her_sc"         "biomass_her_sd_sc"      "biomass_spr_sc"        
#>  [4] "biomass_spr_sd_sc"      "density_cod_rec_sc"     "density_cod_sc"        
#>  [7] "density_fle_rec_sc"     "density_fle_sc"         "density_saduria_rec_sc"
#> [10] "density_saduria_sc"     "depth_rec_sc"           "depth_sc"              
#> [13] "oxy_rec_sc"             "oxy_sc"                 "sigma_E"               
#> [16] "sigma_O"                "temp_rec_sc"            "temp_sc"

# I want the random effects to be dark gray...
pal <- brewer.pal(n = 8, name = "Dark2")
pal2 <- c(pal[1:5], "black", pal[6:8])

# Sort the terms so that random effects are at the top...
mfull_est <- mfull_est %>% 
  mutate(term2 = ifelse(term %in% c("sigma_E", "sigma_O"), 2, 1))

ggplot(mfull_est, aes(reorder(term, term2), estimate, color = group, fill = group, shape = scale)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray40", alpha = 0.5) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
  geom_point(size = 2.5) +
  scale_color_manual(values = pal2, name = "") +
  scale_shape_manual(values = c(15, 19, 21, 17), name = "") +
  scale_fill_manual(values = rep("white", 9), name = "") +
  labs(x = "", y = "Standardized coefficient") +
  scale_x_discrete(breaks = levels(mfull_est$term),
                   labels = c(expression(Herring[rec]),
                              expression(Herring[sub]),
                              expression(Sprat[rec]),
                              expression(Sprat[sub]),
                              expression(Cod[rec]),
                              expression(Cod[haul]),
                              expression(Flounder[rec]),
                              expression(Flounder[haul]),
                              expression(Saduria~entomon[rec]),
                              expression(Saduria~entomon[haul]),
                              expression(Depth[rec]),
                              expression(Depth[haul]),
                              expression(Oxygen[rec]),
                              expression(Oxygen[haul]),
                              expression(sigma[E]),
                              expression(sigma[O]),
                              expression(Temp[rec]),
                              expression(Temp[haul])
                              )) +
  coord_flip() +
  theme(plot.margin = unit(c(0.4, 0.4, 0.4, 0), "cm"),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  guides(color = "none", fill = "none", shape = guide_legend(ncol = 2)) + 
  NULL

  
ggsave("figures/effect_size.pdf", width = 17, height = 17, units = "cm")


# Just a test to see the labels were alright
#ggplot(mfull_est, aes(reorder(term, estimate), estimate)) +
# ggplot(mfull_est, aes(term, estimate)) +
#   geom_hline(yintercept = 0, linetype = 2, color = "red", alpha = 0.5) +
#   geom_point(size = 2) +
#   labs(x = "", y = "Standardized coefficient") +
#   coord_flip()

# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est, "output/mfull_est.csv")

Plot in space

pred_map <- predict(mfull, newdata = pred_grid)

pal <- wes_palette("Zissou1", 21, type = "continuous")

plot_map_raster +
  geom_raster(data = filter(pred_map, depth < 120 & depth > 10), aes(x = X*1000, y = Y*1000, fill = exp(est))) +
  scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
  facet_wrap(~ year, ncol = 5) +
  theme_facet_map() +
  NULL


ggsave("figures/supp/condition/all_years_condition_map_covars.png", width = 6.5, height = 6.5, dpi = 600)

# And a smaller map for selected years
plot_map_raster_labels +
  geom_raster(data = filter(pred_map, year %in% c(1994, 2001, 2008, 2018) & depth < 120 & depth > 10),
              aes(x = X*1000, y = Y*1000, fill = exp(est))) +
  scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") + 
  facet_wrap(~ year, ncol = 2) +
  theme_plot() +
  NULL


ggsave("figures/condition_map_covars.pdf", width = 17, height = 17, units = "cm")



# Plot correlation between random effects and covariates

# Drop NA covariates for the correlation plot
# colnames(pred_map)
# 
# pred_map_corrplot <- pred_map %>%
#   drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
#           density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc, 
#           density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
#           oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc, omega_s, epsilon_st) %>% 
#   dplyr::select(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
#           density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc, 
#           density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
#           oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc, omega_s, epsilon_st)
# 
# corr_re <- round(cor(pred_map_corrplot), 2)
# 
# ggcorrplot(corr_re, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
#   theme_classic() + 
#   theme(axis.text.x = element_text(angle = 90, vjust = 0.3))
# 
# ggsave("figures/supp/condition/correlation_variables_random_effects.png", width = 6.5, height = 6.5, dpi = 600)

Plot spatial and spatiotemporal random effects:

# Plot spatiotemporal random effect
plot_map_raster +
  geom_raster(data = pred_map, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
  scale_fill_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  ggtitle("Spatiotemporal random effects") +
  NULL


ggsave("figures/supp/condition/epsilon_st_map.png", width = 6.5, height = 6.5, dpi = 600)

# Plot spatial random effect
plot_map_raster +
  geom_raster(data = filter(pred_map, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
  scale_fill_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  theme_plot() + 
  ggtitle("Spatial random effects") +
  NULL


ggsave("figures/supp/condition/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)

Calculate the “spatial trend” from the estimates:

# Fit a linear model to each prediction grid of the estimate over time
# https://community.rstudio.com/t/extract-slopes-by-group-broom-dplyr/2751/7
time_slopes_by_year <- pred_map %>%
  mutate(id = paste(X, Y, sep = "_")) %>%
  split(.$id) %>%
  purrr::map(~lm(est ~ year, data = .x)) %>%
  purrr::map_df(broom::tidy, .id = 'id') %>%
  filter(term == 'year')

# Plot the slopes
time_slopes_by_year2 <- time_slopes_by_year %>%
  separate(id, c("X", "Y"), sep = "_") %>%
  mutate(X = as.numeric(X),
         Y = as.numeric(Y))

plot_map_raster +
  geom_raster(data = time_slopes_by_year2, aes(x = X * 1000, y = Y * 1000, fill = estimate)) +
  scale_fill_gradient2(midpoint = 0) +
  facet_wrap(~ term, ncol = 5) +
  theme_plot() +
  ggtitle("Time slopes by each pred grid") +
  NULL

<img src=“condition_model_cf_files/figure-html/calculate”spatial trend”-1.png” width=“672” style=“display: block; margin: auto;” />


ggsave("figures/supp/condition/spatial_trend_condition.png", width = 6.5, height = 6.5, dpi = 600)

Visualize marginal effects

Depth

# Prepare prediction data frame
nd_dep <- data.frame(depth_sc = seq(min(d$depth_sc), max(d$depth_sc), length.out = 100))

nd_dep <- nd_dep %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_dep <- predict(mfull, newdata = nd_dep, se_fit = TRUE, re_form = NA)

mar_dep <- ggplot(p_margin_dep, aes(depth_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Depth[haul]))

Oxygen

# Prepare prediction data frame
nd_oxy <- data.frame(oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100))

nd_oxy <- nd_oxy %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         #oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_oxy <- predict(mfull, newdata = nd_oxy, se_fit = TRUE, re_form = NA)

mar_oxy <- ggplot(p_margin_oxy, aes(oxy_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se)))+
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Oxygen[haul]))

Temperature

# Prepare prediction data frame
nd_tem <- data.frame(temp_sc = seq(min(d$temp_sc), max(d$temp_sc), length.out = 100))

nd_tem <- nd_tem %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         #temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_tem <- predict(mfull, newdata = nd_tem, se_fit = TRUE, re_form = NA)

mar_temp <- ggplot(p_margin_tem, aes(temp_sc, exp(est),
                                     ymin = exp(est - 1.96 * est_se),
                                     ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Temperature[haul]))

Cod

# Prepare prediction data frame
nd_cod <- data.frame(density_cod_sc = seq(min(d$density_cod_sc), max(d$density_cod_sc), length.out = 100))

nd_cod <- nd_cod %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         #density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_cod <- predict(mfull, newdata = nd_cod, se_fit = TRUE, re_form = NA)

mar_cod <- ggplot(p_margin_cod, aes(density_cod_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  theme(aspect.ratio = 1) + 
  labs(x = expression(Cod[haul]))

Sprat

# Prepare prediction data frame
nd_spr <- data.frame(biomass_spr_sd_sc = seq(min(d$biomass_spr_sd_sc), max(d$biomass_spr_sd_sc), length.out = 100))

nd_spr <- nd_spr %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         #biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_spr <- predict(mfull, newdata = nd_spr, se_fit = TRUE, re_form = NA)

mar_spr <- ggplot(p_margin_spr, aes(biomass_spr_sd_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  #coord_cartesian(ylim = c(-4.65, -4.4)) + 
  theme(aspect.ratio = 1) + 
  labs(x = expression(Sprat[sub]))

Saduria

# Prepare prediction data frame
nd_sad <- data.frame(density_saduria_sc = seq(min(d$density_saduria_sc), max(d$density_saduria_sc), length.out = 100))

nd_sad <- nd_sad %>%
  mutate(year = 2003L,
         year_f = as.factor(2003),
         ln_length_cm = 0,
         biomass_her_sc = 0,
         biomass_her_sd_sc = 0,
         biomass_spr_sc = 0,
         biomass_spr_sd_sc = 0,
         density_cod_sc = 0,
         density_cod_rec_sc = 0,
         density_fle_sc = 0,
         density_fle_rec_sc = 0,
         #density_saduria_sc = 0,
         density_saduria_rec_sc = 0,
         depth_sc = 0,
         depth_rec_sc = 0,
         oxy_sc = 0,
         oxy_rec_sc = 0,
         temp_sc = 0,
         temp_rec_sc = 0
         )

# Predict from full model (AIC-selected)
p_margin_sad <- predict(mfull, newdata = nd_sad, se_fit = TRUE, re_form = NA)

mar_sad <- ggplot(p_margin_sad, aes(density_saduria_sc, exp(est),
                                    ymin = exp(est - 1.96 * est_se),
                                    ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  #coord_cartesian(ylim = c(-4.75, -4.4)) + 
  theme(aspect.ratio = 1) + 
  labs(x = expression(Saduria~entomon[haul]))

Plot together

# mar_dep + mar_oxy + mar_temp + mar_cod + mar_spr + mar_sad
#   plot_annotation(tag_levels = 'A') 
# 
# ggsave("figures/supp/condition/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)

p_margin_dep2 <- p_margin_dep %>%
  mutate(variable = "Depth (haul)") %>% 
  dplyr::select(est, est_se, depth_sc, variable) %>% 
  rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_oxy2 <- p_margin_oxy %>%
  mutate(variable = "Oxygen (haul)") %>% 
  dplyr::select(est, est_se, oxy_sc, variable) %>% 
  rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_tem2 <- p_margin_tem %>%
  mutate(variable = "Temperature (haul)") %>% 
  dplyr::select(est, est_se, temp_sc, variable) %>% 
  rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_cod2 <- p_margin_cod %>%
  mutate(variable = "Cod (haul)") %>% 
  dplyr::select(est, est_se, density_cod_sc, variable) %>% 
  rename(value = density_cod_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_spr2 <- p_margin_spr %>%
  mutate(variable = "Sprat (sub)") %>% 
  dplyr::select(est, est_se, biomass_spr_sd_sc, variable) %>% 
  rename(value = biomass_spr_sd_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_sad2 <- p_margin_sad %>%
  mutate(variable = "S. entomon (haul)") %>%
  dplyr::select(est, est_se, density_saduria_sc, variable) %>% 
  rename(value = density_saduria_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

margs <- bind_rows(p_margin_dep2, p_margin_oxy2, p_margin_tem2, 
                   p_margin_cod2, p_margin_spr2, p_margin_sad2)

ggplot(margs, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
                  ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  facet_wrap(~variable, scales = "free_x") +
  labs(y = "Le Cren's condidition factor",
       x = "Scaled value") +
  theme_plot() +
  theme(axis.text.x = element_text(angle = 0)) +
  NULL


ggsave("figures/supp/condition/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)